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Abstract 

We present spectral methods developed in our group to solve three-dimensional 
partial differential equations. The emphasis is put on equations arising from astro- 
physical problems in the framework of general relativity. 



1 Introduction 

Astrophysical problems involve many field of physics: hydrodynamics, mag- 
neto hydrodynamics, chemistry, kinetical, nuclear and statistical physics, grav- 
itation... Physicists and astrophysicists have developed well suited techniques 
to solve numerically mathematical problems in each of these specific fields. 
Moreover they have to handle problems in which different branches of physics 
interact and therefore compatibility of different numerical techniques is re- 
quired. 

The following examples illustrate the above statement. Chaotic chemical re- 
actions can occur in the interstellar medium [34]. Moreover the interstellar 
medium is turbulent and mildly supersonic. Therefore the numerical mod- 
elizations requires numerical methods able to describe chaos and compressible 
turbulence. Astrophysicists involved in numerical work are confronted with 
problems that require the integration of systems of quasi-linear partial dif- 
ferential equations (PDEs) (i.e. PDEs which are linear in the highest order 
derivatives). In general these systems of equations contain all kind of PDEs 
i.e. hyperbolic, parabolic and elliptic PDEs. The study of stellar oscillations 
is another good example. The system of PDEs that must be integrated is the 
Navier-Stokes equation for the fluid motion coupled to the continuity equation 
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for the mass conservation, to the scalar Poisson equation for the gravitational 
field and to the heat equation (since the plasma of the star has a finite heat 
conductivity) . 

PDEs systems for special relativistic astrophysics present the same structure 
than the Newtonian ones but special relativistic hydrodynamic equations are 
more complicated because of the presence of high Lorentz factors which make 
life less easy. Relativistic jets in active galactic nuclei are the realm of spe- 
cial relativistic [magneto-] hydrodynamics codes. Things change drastically for 
compact astrophysical objects such as neutron stars and black holes. For these 
objects the Newtonian theory of the gravitation is no longer appropriate and 
general relativity must be used. 

Differential geometry is the realm of GR. New geometrical concepts such as 
topology, globality, isometry, Lie derivative, Killing vector, curvature, etc... 
play an essential role. Of course all these concepts are present in the frame- 
work of special relativity, but in a trivial and underlying way. Therefore, as- 
trophysicists dealing with special relativity, as a modern Monsieur Jourdain 
who "fait de la prose sans le savoir" (he speaks in prose without knowing it), 
use these concepts without realizing it. 

In general relativity, the gravitational potential is no longer scalar. The metric 
tensor g a p, (a, (5 = 0,1,2,3) replaces it, its 10 independent components playing 
the role of gravitational potentials. The gravitational field equations (Einstein 
equations) can be written in a beautiful very compact form: 

Ra(3 ~ ^Rgaf3 = T a/3 , (1) 

where R a p = Rp a and R are respectively the Ricci tensor and the curvature 
scalar R = R a a associated with the metric g a p, and where T a/3 is the energy- 
momentum tensor. The equations of motion can be written in a compact 
way too : V ' a T a p = 0. If an electromagnetic field is present, the Maxwell 
equations must be added to the above equations : V a F a p = Jp where F a p is 
the electromagnetic tensor and Jp is the quadri-electric current. The simplicity 
of these equations is only apparent. In fact, when explicited, they contain a few 
thousand of terms. Finding a solution to a complete 3-D physical problem in 
the framework of general relativity is the Holly Grail of Numerical Relativity. 

Einstein equations form a system of 10 second order quasi-linear PDEs. More- 
over, the 4 identities (Bianchi identities) 

V a (R% - l -Rg a p) = (2) 
introduce 4 degrees of freedom for the solution. This is equivalent to the free- 
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dom of gauge choice in electromagnetism. The choice of the gauge in GR 
defines the character (i.e. ellipticity, parabolicity or hyperbolicity) of the equa- 
tions. Actually, the Einstein's equations are degenerated in the sense that their 
character depends on the gauge choice. For example, these equations are all 
hyperbolic when the so-called harmonic gauge [1] is chosen and, on the con- 
trary, the system can be reduced to 5 elliptic equations coupled to 5 hyperbolic 
equations if the so-called radiation gauge [53] is chosen. 

The first problem which arises naturally when one wants to solve the Einstein 
equation for a particular astrophysical problem is the determination of the 
most appropriate gauge. There is no definitive answer to this question and the 
debate is still open. One has to consider essentially two criteria. The first one 
depends on the physics of the problem: depending on the coordinate system 
the space-time slicing induced by the gauge may or may not be able to describe 
and cover the physical event that one wants to study. The second one depends 
on the numerical technique employed: depending on the gauge one may or 
may not be able to find a numerical solution to the equations (e.g. ability to 
solve coupled non-linear elliptic equations). It is tempting to choose a gauge 
which seems more appropriate to the numerical technique used. However, we 
don't think that this is a good approach to the problem. 

We hope to have given to the reader a vague idea of the kind of problems 
met by astrophysicists and of the mathematical difficulties that have to be 
overcome. We will show in this paper how spectral methods may be used to 
successfully overcome these difficulties. Spectral methods (SM) will be intro- 
duced in details in Sect. 2. For now, let us anticipate two of the most important 
advantages of SM. The first one is the economy of the number of degrees of 
freedom (NDF). For a given numerical accuracy, less degrees of freedom are 
required than in the case of finite difference techniques. More precisely, as a 
rule of thumb, we can say that with SM, the number of grid points is reduced 
by a factor of five per space dimension. This means that, in the case of 3-D 
problems, the number of grid points will be decreased by a factor of 5 3 = 125. 
Now, if one consider a dynamical problem where the integration time step 
scales at least as N, where N is the NDF, this advantage becomes obvious. 
The second advantage consists in the possibility of using a coordinate sys- 
tem well adapted to the geometry of the problem and to handle exactly the 
pseudo-singularities potentially present in the chosen coordinate system. For 
instance, thanks to SM it has been possible to study 3-D turbulent motion 
of an incompressible fluid with Reynolds number of about 1000 [55] and to 
describe the disruption of a star in the tidal field of a black hole in less than 
3 hours of C.P.U. time on a workstation [39] when the same problem, using 
a finite difference method, has required 30 hours of C.P.U. time on a Cray-2 
[32]. Moreover, thanks to the high precision achievable with SM the stability 
and collapse of a N.S. close to its critical mass was successfully studied in the 
GR frame [22], [26]. 
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Elliptic and parabolic PDEs are the realm of SM. Nevertheless good results 
can be obtained for hyperbolic equations too if no discontinuities appear in 
the solution [16], [39]. When a discontinuity is present, some kind of natural 
or artificial viscosity must be added in order to obtain a sufficiently smooth 
solution. Interesting results were obtained with this technique in describing 
a mildly supersonic turbulence occurring in star formation [15,48,54,49], or, 
by our group, in modelling a 3-D stellar collapse [37]. Other numerical meth- 
ods (especially those based on Riemann solvers technique [40]) have been 
developed in order to handle shocks. Successful results were obtained in mod- 
elizing relativistic jets (see [31] and literature quoted there). In the case of 
relativistic problems where systems of elliptical coupled to hyperbolic PDEs 
must be solved, an interesting strategy consists in using SM for the elliptical 
equations and Riemann solver technique for the hyperbolic equations where 
discontinuities may occur. Interesting results have already been obtained in 
1-D modellisation of the collapse of a core of a supernova in the framework 
of a tensor-scalar theory of gravitation [42], [43]. This strategy looks very 
promising too for 3-D problems. 

Our aim in the present paper is to explain how SM work and how they can be 
implemented in the easiest way. The mathematical theory of SM will be left to 
specialized literature (e.g. [21], [18]). The reader does not need to know astro- 
physics but almost all the examples presented as illustration are taken from 
the modelization of astrophysical phenomena. Therefore some short explana- 
tion of the phenomenon that has motivated the modelization will be given. 
On the contrary some elementary notions of differential geometry would be 
worth for understanding the mathematical problems. 

The plan of this paper is as follows. We give first a short introduction to 
spectral methods (Sect. 2) in which the essential notions are presented. Then 
we explicit the algebraic properties of representation of some relevant dif- 
ferential operators (Sect. 3). Basic methods for solving PDEs are presented 
in Sect. 4, especially the time scheme and the treatment of boundary con- 
ditions. This is illustrated by selected one-dimensional examples (Sect. 5). 
Three-dimensional (3-D) spectral methods with spherical-like coordinates are 
introduced in Sect. 6. Scalar and vectorial equations, such as Poisson or tele- 
graph equations, are considered in Sects. 7 and 8 respectively. After a brief 
discussion of the structure of Einstein equations, the paper ends by the pre- 
sentation of various astrophysical results (Sect. 9). 



2 Introduction to spectral methods 
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2.1 A simple example based on a Fourier expansion 



In this section we present the main ideas underlying spectral methods (SM). 
We shall give only the basic notions, living the theory to specialized and very 
good existing text books, among which we recommend the monographies by 
Gottlieb & Orzsag [21] and by Canuto, Hussaini Quarteroni & Zang [18]. 

Let us consider the quasi-linear partial differential equation (viscous Burger 
equation): 

du d 2 u , du . 



where u is a function of the two variables t and x and A is some parameter. 

Consider first the simple case A = (linear heat equation). Let us assume that 
the value of u(t, x) is known for all x at t — and that u is a periodic function. 
The basic idea underlying SM consists to transform the PDE in a system of 
ordinary differential equations by means of a expansion of the solution onto a 
series on a complete basis. Since u is assumed to be periodic, it is natural to 
use a Fourier expansion: 

oo 

u(x, t) — [dk(t) cos (2nkx) + bk(t) sin {2-Kkx)\ . (4) 

fc=0 



Equation. (3), with A = 0, can then be rewritten as 



In this way, finding a solution of the PDE (3) turns out to be equivalent to 
solve the infinite system of ordinary differential equations (5), where the initial 
values of the coefficients and bk are given by the Fourier expansion of the 
function u at the time t — 0: 

i i 
Gfc(0) = J u(x, 0) cos(27r/c:r) dx, b k (0) — J u(x, 0) sin(27rA;x) dx . (6) 





From the numerical point of view, the series (4) has to be truncated. Let be 
N/2 — 1 the highest term of the series (N even). The integrals (6) must be 
evaluated numerically in the most accurate way: if u(x, 0) does not contain 
spatial frequencies higher than N/2 — 1, then its Fourier coefficients a k and 
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bk must be computed exactly within the roundoff errors. In other words we 
require that the numerical integrals 




J sin(27rA;x) cos(2irjx)dx 
o 



(7) 



are computed exactly. The minimum number of grid points at which the 
functions u(x, 0) must be sampled is N. The sampling points that fulfil the 
above requirements are called the Gauss-Lobatto points or collocation points: 
Xj = j/N. These points are equally spacedQ. 

Note that there exists an isomorphism between the coefficients a&(t), and 
the values u(xj,t) of the function at the collocation points. For this reason, 
N is called the number of degrees of freedom (NDF), denomination that we 
prefer to the more usual number of grid points. Finally it is worth to note that 
the uniform repartition of the Gauss-Lobatto points should be considered as 
an accident. If one changes the basis of expansion, uniform spacing does not 
a priori correspond to a Gauss-Lobatto sampling (see the next section). 

Consider now the non linear problem A ^ 0. We can proceed in two different 
ways to treat the non-linear term. The first one (spectral method in the strict 
sense) consists in computing the first space derivative of u in the Fourier space 
(or coefficient space) and to perform a convolution with Uk- From a numerical 
point of view, these terms, computed at the time P, act as a source and allow 
to compute the coefficient Uk at the time P +l no matter the temporal scheme 
used. Within spectral method, all the evolution of the equation is computed 
in the Fourier space. Backing in the configuration space (at arbitrary values 

1 The properties of the Gauss-Lobatto sampling for the discrete Fourier decompo- 
sition can be recovered easily by bearing in mind that J2n=d 1 n = 0- ~ 1 N )/(^ ~ 
<?)> \q\ < b from which the following identity holds (q = exp(27ri(k — j)/N)): 



This means that the computation of the integrals (7) by means of the trapeze 
method is exact (within the round-off errors). In the theory of signal, the Gauss- 
Lobatto points correspond to the Shanonn sampling, and N/2 — 1 is related to the 
Shanonn frequency. 



i N-l 
N / exp(2i7r(A; — j)x)dx = exp(2nr(£; 



3)Xn 



N-l 

/J exp(27ri(/c 

n=0 



j)n/N) = N5 kj 



(8) 
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of x) is performed only for practical convenience, for example to visualize 
the solution at a given time. Consequently, the configuration space plays a 
secondary role. 

The second way to proceed (pseudo-spectral method) consists in computing the 
derivative of u in the Fourier space, to come back in the configuration space 
by an inverse Fourier transform, to multiply du/dx by u in the configuration 
space and to come back again in the Fourier space. With pseudo- spectral 
method, the Fourier and configuration spaces behave equally. Pseudo-spectral 
method is generally used in spite of the apparent athletic work of dancing from 
one space to the other. The reason is that performing a convolution requires 
a number of operations proportional to iV 2 . On the contrary, thanks to the 
fastQ Fourier algorithm, the number of operations needed by pseudo-spectral 
method is proportional to iVlogiV. Moreover the pure spectral method can 
not handle easily non-linearities more general than the quadratic one present 
in (3). In all this paper, we shall use (incorrectly) for shortness the term of 
spectral method instead of the more appropriate pseudo-spectral method. 

The above example contains all the philosophy of SM. A PDE, or more gen- 
erally a system of PDEs, is transformed in an equivalent system of ordinary 
differential equations. The set of trigonometric functions used for the expan- 
sion, has been chosen because it fulfils automatically the boundary conditions 
and because it exists a fast transform algorithm. The space derivatives are 
computed in the coefficients space and the computation uses the value of the 
function on all collocation points. Quadratic terms are computed by coming 
back to the configuration space. The most important consequence is that if 
the solution is C°° then the C 2 error between the numerical solution and the 
analytical one scales as exp (-N) and as aMp+i) f or a solution C*Q. An error 
that scales as exp(— N) is called evanescent. 

Finally the computation of the required quantities is performed in one of the 
two spaces, namely the configuration space or its dual (Fourier space) , in order 
to save computation time and to obtain high accuracy. 

2 A 1-D numerical algorithm is called fast if the number of arithmetic operations 
scales with the number of the degrees of freedom N no faster than N log N. Fast 
Fourier Transform (FFT) and Fast Chebyshev transform algorithms (FCT) are fast 
algorithms. 

3 Following the terminology of the functional analysis, the error L v of some numer- 
ical solution u num (x,t) is defined by 




0) 



where u cx (x,t) is the exact solution. 
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N = 4 



JV = 8 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 



x x 

Fig. 1. Graph of the function u(x) = cos 3 (7rx/2) — (x + l) 3 /8 (solid line) and its 
Chebyshev interpolant (dashed line) for N = 5 (left) and N = 9 (right). For N = 9, 
the two curves cannot be distinguished graphically. The circles denote the values 
at the collocation points x n . Note that spurious oscillations between the collocation 
points are not present. 



iV = 4 N = 8 




-1 -0.5 0.5 1 -i -0.5 0.5 1 



x x 

Fig. 2. Graph of the derivative of the function represented in Fig. 1 Eq. (12) (solid 
line) and its Chebyshev approximation (dashed line), for N = 5 (left) and N = 9 
(right). The circles denote the values at the collocation points x n . 



8 



2.2 More general expansions 



If the initial conditions or the boundary conditions are not periodic, Fourier 
expansion is no more adapted to treat the problem because of the presence of 
a Gibbs phenomenon at the boundaries of the interval 0. A new base of func- 
tions must be used. Expansions in Legendre or Chebyshev polynomials are the 
most common ones. The reader will find in [21] and [18] a discussion on the 
advantage and drawbacks of the different expansions. Legendre polynomials 
are often used for theoretical work [36]. The main advantage of Chebyshev 
polynomials is that there exists a fast algorithm to perform Chebyshev expan- 
sion (its C.P.U. time is just a little bit longer than the corresponding Fourier 
transform time). For this reason, Chebyshev expansion is widely used in this 
paper. 

Chebyshev polynomials T n (x) are defined as the family of polynomials which 
obeys to the the following orthogonality relations: 

(1 + <W f T n (x)T m (x) 

/ — 7= = ax = b nm , (10) 

7T J v 1 — x 2 



where 8 nm is the Kronecker symbol. By introducing the new variable 9 = 
arccos(x) (0 < 9 < n), the Chebyshev polynomials read 

T n (x(6)) = cos(n0) . (11) 



The associated Gauss-Lobatto points Xj are uniformly spaced with respect to 
the variable 9: Xj = cos(ir(j — 1)/(N — 1)), < j < N — 1 where N is the 
number of degrees of freedom. It easy to see that the collocation points are 
more dense close to the boundaries of the interval [— 1,+1]. Thanks to the 
relation (11), Chebyshev polynomials have close links to the trigonometric 
functions used in the Fourier expansion. 

The following example illustrates the accuracy of SM. Let us consider the 
function 

u(x) = cos 3 (^x^j -i(s + l) 3 x G[-l,l]. (12) 



4 The Gibbs phenomenon occurs when an expansion, for example the Fourier ex- 
pansion, of a discontinuous function is performed. In this case the C 2 error scales 
with the degrees of freedom N as 1/N. On the contrary the pointwise error C°° 
tends to a finite value when N — > oo (see e.g. ref. [21]) 
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Table 1 

Comparison between a Chebyshev spectral method (SM) and a second-order finite 
difference method (FDM) in the computation of the first and second derivatives of 
the function u{x) = cos 3 (|x) — (x + l) 3 /8. The error is defined as the highest dis- 
crepancy with the exact value at the grid points. 16 digits are used in the numerical 
computation. Note that to obtain an accuracy on d 2 u/dx 2 of the order 10~ 4 , 13 
points are sufficient for the SM, whereas the FDM requires about 500. Note also 
that for NDF > 129, the round-off errors make the SM less efficient. 

Note that u(x) is not a polynomial, so that its expansion over the Chebyshev 
polynomials a priori requires an infinite number of coefficients. The term (x + 
l) 3 has been added to the cosine term in order not to have a periodic function. 
The function u(x) is represented in Fig. 1, as well as its Chebyshev interpolant 
Y,n=o u n T n (x) for iV = 5 and N = 9. Note that for a number as small as 
N = 9, the Chebyshev expansion cannot be graphically distinguished from 
u. The derivative of u(x) is numerically computed by means of the method 
described in Sect. 3.1 and is compared with its analytical value in Fig. 2. 

The comparison with a second-order finite difference method is performed in 
Table 1 for the computation of the first and second derivatives of u. Note that 
for N = 33, the SM reaches a relative accuracy of 1.5 x 1CT 13 and 3 x 1CT 11 
in the computation of du/dx and d 2 u/dx 2 respectively. 
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3 Algebraic properties and spectral representation of some ele- 
mentary operators 

We will discuss in this section some properties and implementation of elemen- 
tary linear differential operators when using Chebyshev polynomial expan- 
sions. 



3.1 First and second order derivatives. 



The expansion coefficient bi of df jdx are given by 6j = D^aj where are the 
coefficients of / and, by definition of d/dx, 

1 

Dij = (2 - <Sj) J T^T'jix^l - x 2 )- l ' 2 dx . (13) 
-l 



The matrix is computed from the relation 

n + l 

T' n+1 (x) = 2(n + l)T n (x) + —-T' n _ 1 (x) , (14) 

71/ X 



which can be derived easily from Eq. (11). For a NDF = 7, the result is: 
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(15) 



The disposition of the non vanishing elements of is due to the property 
of the operator d/dx which transforms the parity of the functions at which is 
applied. The main diagonal and the last line of vanish because the operator 
d/dx transforms a polynomial of degree n into a polynomial of degree n — 1. 
These two important properties will be used below. 
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The representation of the operator d 2 /dx 2 can be obtained straightforwardly 
by taking the square of the matrix (15). For a NDF = 7, it reads 
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The disposition of the zeros is different, because the operator d 2 /dx 2 preserves 
the parity. The last two lines of the matrix vanish because this operator de- 
creases the degree of a polynomial by two units. Other operators that preserve 
the number of coefficients are 

x l- (ax + bx2) & (17) 

These operators play an important role in solving PDEs. 

The computation of the derivative of a function / obtained by a matrix mul- 
tiplication needs a number of multiplication oc iV 2 . Actually, the particular 
structure of the matrix (15) suggests not to perform the full matrix multipli- 
cation but instead to use the recursion formula 

f> n = 2(n + l)f n+1 + f> n+2 n<N-2, (18) 

A similar recursion law exists for the operator d 2 /dx 2 . In a such a way the 
number of operations becomes oc N. Hereafter we shall call fast an algorithm 
requiring a number of elementary operations oc iV or oc N In N where N is the 
NDF. 



3. 2 Primitive of a function 



A primitive F(x) of a function f(x) can be obtained by inversion of the oper- 
ator d/dx, i.e. by solving the algebraic system of N equations 

N-l 

^ ' Dmnbn &m j (19) 
n=0 

where b m are the expansion coefficients of F . The determinant of D nm vanishes 
because the primitive of a function is defined within some additive constant. 
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We can add an algebraic condition, for instance F(—l) = 0, in order to elimi- 
nate this degeneracy. One (certainly not unique !) way to proceed is to seek for 
a primitive that has some property in the coefficient space, for instance that 
its first coefficient b vanishes. That is equivalent to find a particular solution 
F p of the differential equation dF/dx = f. Once this solution is found, we can 
add to it a homogeneous function Fh in such a way that the final solution 
F = F p + Fh fulfils the required conditions. In this example Fh is constant. If 
we look for a particular solution having b = 0, the new system is 



1 






















' b 










2 



















h 




a 








4 
















b 2 




di 











6 













h 




d 2 














8 














d 3 

















10 







h 




d± 




















12 








a~5 



(20) 



where dj is the new RHS vector corresponding to the linear combination 

Dmn = Ann - — ^ An+2,n < m < iV — 4 , (21) 

' m 



which transforms the matrix D mn into a diagonal one. It follows that the com- 
putation of the primitive of a function can be performed with a fast algorithm. 
This technique is quite general and is similar to that used when treating more 
complicated problems (cf. Sect. 4.5). 



3. 3 Reduction to band matrices of more general operators 



The matrix representation F) 2 mn of the operator d 2 /dx 2 can also be transformed 
to a band matrix with some simple linear combination of the lines. Applying 
first the linear combination 

D 2 = U+OoMAnn ^n+2,n Q < m < iV - 3 (22) 

m + 1 ~ v ' 



followed by 

D 2 mn = D 2 mn -D 2 m+2n 0<m<N-5 (23) 
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transforms the matrix into 
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D 2 = 

nm 
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More generally, the matrix representation of the differential operator 

d d 2 
L = I + (a + aix) — + (&o + hx + b 2 x 2 ) — , (25) 

where I is the identity and ao, on, bo, b\ and 62 are numerical constants, can 
be transformed to a penta-diagonal matrix via the linear combinations (22) 
and (23). Therefore there exists a fast algorithm to inverse the above operator. 
Note that the operator I is responsible for the five diagonals of the matrix. 

Finally, the multiplication and division of a function by X = a + bx can be per- 
formed in the coefficient space. The matrix representation of this operation is 
a bi-diagonal matrix X^. Therefore multiplication and division are performed 
with a fast algorithm. The advantages in operating in this way are the follow- 
ing ones. Firstly, it is easier to handle the division by X when X vanishes for 
a given value of x. Moreover, the roundoff error is less important than when 
performing the multiplication in the configuration space. Secondly, computing 
quantities like (x + l) 1 , (I integer > 1) in the configuration space introduces 
high frequency terms. In fact, for large values of /, the function (1 + x) 1 is 
numerically zero and therefore not continuous near x = — 1. Such functions 
play an important role in spherical coordinates, because they are the radial 
part of the spherical harmonics of order /. Thirdly, it avoids to come back to 
the configuration space to perform the multiplication by x when the function 
is known in the coefficient space. 



4 Solution of elementary P.D.E. in Cartesian coordinates 



4-1 Temporal scheme 



Consider the 1-D second order hyperbolic PDE (heat equation with advection) 
(90 (90 (9^0 

lk +v{x) ix~~^ x) ^ = S{x ' t) ^[-1,+1], t>0. (26) 
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If /i(x) = and v(x) — vq — const, the general solution is Q(x, t) = f(x — vot) 
i.e. a wave propagating from the left to the right if vq > and from the right 
to the left in the opposite case. The time evolution of 6 can be obtained by 
computing the time derivative with a finite difference method: t is discretized 
at given uniformly spaced instants tj, tj+i = tj + At and a second order scheme 
(for example an explicit Crank-Nicolson scheme) approximates, Eq. (26) by 

where dQi +1 /dx at the RHS is computed by extrapolating its values from the 
past. This scheme is highly unstable, no matter the choice of the time-step 
At. This instability is due to the absence of boundary conditions. The value 
of at x = — 1 if Vq > or at x = 1 in the opposite case, must be given. In 
other words, we have to define what goes into the interval [—1, +1]. 

Before to explain which boundary condition makes the problem well posed 
(in the sense given by Hadamard, see e.g. [50] p. 62) and how it can be im- 
plemented, it is worth to discuss the Courant conditions that garanties the 
numerical stability. 

In the case of Chebyshev polynomial expansion, the numerical stability is 
achieved if At oc 1/N 2 . This condition (Courant condition) is much more se- 
vere than in the case of finite difference method with uniform grid for which 
At oc 1/N. The situation is much worse if a spatial second order operator 
is present ^ 0). In this case At must decrease as 1/N 4 . These severe 

Courant condition is due to the fact that the density of the Chebyshev col- 
location points behave like 1/N 2 near the boundaries x = — 1 and x — +1 
(see e.g. Fig. 1 or Fig. 3). However, if the coefficients of the space derivatives 
vanish at the boundaries, the Courant condition becomes At oc 1/N in the 
case of a first order spatial differential operator and At oc 1/N 2 in the case of 
a second order spatial differential operator. Note that the dense sampling at 
the boundaries can also be an advantage: it is well suited to resolve bound- 
ary layers if they are present in the solution. We shall give an example below 
(Sect. 5.2). 

Implicit (or semi-implicit) schemes allow to overcome this difficulty Implic- 
itation of the Crank-Nicolson scheme (27) can be obtained by inverting the 
operator Id + Atd/dx. The generalization to the case fj,(x) — /io — const. ^ 
is straightforward. With an implicit scheme the absolute convergence (for any 
value of At) is obtained. 

The leapfrog scheme, which is very widely used with Fourier expansion, turns 
out to be unstable with Chebyshev expansion. The reader will find an ex- 
haustive discussion on this subject in [21] pp. 103-115. In the present arti- 
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cle, we shall use systematically the Crank-Nicolson scheme. Its accuracy be- 
haves like At 3 per time-step when all the quantities are computed at the time 
P+ 1 / 2 = P + At/2. 

In a more general case, when the coefficients of the spatial derivatives de- 
pend on x, the implicitation consists to invert an operator of the kind O = 
Id + b(x,t)D where Id and D are respectively the identity and a differential 
operator and where b(x, t) is some function. The matrix representation of op- 
erator O is in general a full matrix which cannot be reduced to a band matrix 
by means of simple operations. Consequently, the inversion of O operator re- 
quires a number of operations oc iV 3 . In order to save CPU time and to build 
a fast algorithm, the following semi-implicit method is highly recommended. 

Consider for simplicity Eq. (26) with v(x) = 0. If the Crank-Nicolson scheme 
is used, the corresponding PDE reads 

i f) 2 (~)i +1 1 d 2 P)i 

- 2 At M*) = 05 + 2 At + S W +1/2 ■ (28) 



The above time discretization can be reformulated in the equivalent form 



i d 2 G j+1 1 d 2 <r) j 

e*. =e , + _* Ml) _ 

i d 2 <~)i +1 
--At(^ - I*{x))-qj- + S(xy +1 / 2 , (29) 

where JJj IS the maximum value of fi(x) on the interval [—1,-1-1]. The 
RHS of Eq. (29) becomes the source of the new equation. The term (fi max — 
^i{x))d 2 Q^ +1 / dx 2 is computed by means of a second order extrapolation from 
its value at the items P~ 2 ,P~ X ,P '. In this way, we reduce the problem to the 
case of constant coefficients. 

It is obvious that the new temporal scheme is of second order too. However, its 
precision depends on how close /i(x) is to \i max . Moreover, this method cannot 
be applied straightforwardly to the case where ji{x) vanishes at the boundary 
of the interval. 

To improve the method, we can introduce a polynomial Pi{x) = do + a\x + 
a 2 x 2 instead of fi max which behaves as /i(x) at the boundaries and such that 
Pzix) ~ M^) — 0- A s mentioned in Sect. 3.3, the matrix representation of 
the operator Id + P2{x)d 2 / dx 2 can be reduced to a penta-diagonal matrix 
by means of the linear combination (22)- (23). The generalization to the case 
v(x) 7^ is obvious. 
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Implicit or semi-implicit temporal scheme in multi-dimensional problems can 
be easily implemented with the Alternating Direction Method (cf. [51]). We 
have tested this method for 2-D problems and did not found any difficulty 
We conjecture that this method can be used for 3-D problems too. 

4-2 Boundary conditions and well posed problems 

Because of its generality we will use equation (26) to discuss in detail the 
numerical problems tied to boundary conditions and its solution. 

The problem is well posed in the following cases. 

• If fi(x) = 0, the equation reduces to a first order PDE. In this case, if 
v(x) 7^ for x G [—1, +1] then one BC at x — — 1 ( v(x) > 0) or at x = +1 
(v(x) < 0) must be imposed. 

• If n(x) = and v(x) > for x G [—1, +1] and v(— 1) = then no BC can 
be imposed. The same holds if v(x) < for x G [—1, +1] and v(+l) = 0. 

• If fi(— 1) = yu(+l) = v(x) = but /jl(x) > for x g] — 1, +1[, the rules given 
for the first order PDE hold. 

• If n(x) and v(x) vanish on the boundaries, only one BC is required (as in 
the previous case). 

• If fj,(x) > everywhere, two BC are required. 

For a second order PDE, boundary conditions are quite general and can have 
the form 

a(f)|? + /?(f)e = 7(*) , (30) 

but other non local BC (dual BC) can be chosen (see below). An example of a 
not well posed problem is to impose two BC (the value of 6 and its derivative) 
at one boundary of the interval. 

4-3 Boundary layer and dual boundary conditions 

Consider Eq. (26) with ji and v constant: fx(x) = /i , v(x) — v > 0. If 
Ho — > 0, the equation degenerates in a first order equation. We can not impose 
two BC in the case of a first order equation. Physically, when /x — > 0, a 
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boundary layer is formed at x — +1 and the numerical algorithm has been 
able to describe the strong variation of the solution near the boundary. We can 
form an adimensional parameter N p (the Peclet number) which characterizes 
the thickness of the layer : N p = v L/fi where L is the typical dimension 
of the system. The thickness A of the layer is A = L/N p . In this example, 
L = 2 and N p = 2vq/hq. Thanks to the accumulation of collocation points at 
the boundary, Chebyshev expansion is well suited to treat this problem. The 
layer is well resolved and the maximum value of N p which garanties numerical 
stability scales as iV 7 / 4 . For instance with N = 33 collocation points N p can 
be as high as 200 (see e.g. [21] p. 140). 

In physical applications, it can happen that a fine resolution of the boundary 
layer is not required. The natural question arises how modify the BC in such 
a way that the new solution differs appreciably from the exact one only in the 
layer. It is clear that the boundary layer is created by the BC that not would 
exist in the case /i = 0. Unfortunately, numerical instabilities appear if only 
one BC is imposed when /i 7^ 0. To overcome this difficulty, we need to find a 
criterium which leads to a value of at the boundary where the layer forms 
(x = +1) such that the numerical scheme is stable if /xo = too. In other 
words, we need a redundant BC. 

Such a BC is more general than the usual one. We will call it a "non local" 
boundary condition because this condition is imposed in the coefficients space. 
There is a wide class of criteria to find a magic value for 6(+l). One criterium 
consists in finding 6(+l) such that, at each time step, Q(x) is as smooth as 
possible. For instance, one can minimize the norm of d 2 Q(x)/dx 2 i.e. 



Another method consists to cancel the last coefficient of the expansion of 
or to minimize the sum of the square of the last J coefficients. All these exotic 
BC garanties the stability of the numerical scheme for any value of /io > 
and gives very similar results. Note that the last two criteria introduce an 
evanescent error when /i(x) = 0. 

An heuristic justification of this procedure can be easily understood by noting 
that the Peclet number parametrizes the antagonism between the advection 
term vq > 0, which carries the information from the left to the right, and 
the diffusion term /i which diffuses back the information. For high Peclet 
numbers, the diffusion term influences only the boundary layer. Therefore, 
the solution outside of the layer is almost not influenced by the BC. 

The above exotic BC are useful in other physical situations. Consider a region 




(31) 
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R of the interior of a star. Let us suppose that we are interested in phenomena 
described by a PDE of the kind of Eq. (26) only in this region. We need to 
know which BC to impose at the boundary of R. This BC cannot be known 
without solving the equation in all the star. We know only that the solution 
is smooth in the interior of the star. Therefore, we can impose a value of the 
solution which makes the solution as smooth as possible. Roughly speaking, 
the underlying philosophy is to introduce the minimum of information when 
the problem is not completely determined. The reader will recognize the same 
philosophy which underlies the maximum entropy criterium used in the theory 
of signal. 

4-4 Implementation of the boundary conditions 

There are various ways to implement the BC within spectral methods. We 
present here three different methods. 

4-4-1 Galerkin method 

The Galerkin method consists in choosing the set of functions used for the 
expansion a set of function satisfying the required BC. The example showed in 
Sect. 2.1 solves a PDE with periodical BC by means of Fourier expansion. The 
Galerkin method will be widely used to handle regularity conditions (e.g. with 
spherical coordinates where coordinate singularities appear). It is sometimes 
easy to find a set of functions satisfying the desired BC by taking a linear 
combination of Chebyshev polynomials. For instance, if the solution of the 
PDE must vanish at x — —1, one can expand the quantities onto the new 
basis of function <£>„ : 

$ n (x) =T n (x)+T n+1 (x) . (32) 

Note that the matrix of the corresponding linear algebraic system of equa- 
tions can be reduced to a band matrix via the linear combinations (22), (23). 
However, the number of diagonal is larger: the 5 diagonals for Eq.(26) become 
7 when the above Galerkin basis is used. 



4-4-2 Lanczos method 

The Lanczos method [33] consists to replace the source S 3 n +1/2 (in the coefficient 
space) of Eq. (26) by 

S n = St 1 ^ + MiV-l,n + MiV,n , (33) 
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where 5 n ^ m is the Kronecker symbol and to compute, at each time step, the 
coefficients b\ and b 2 i n such a way that BC are fulfilled (cf. Sect. 4.5). 

4-4-3 Dual Lanczos method 

A dual Lanczos method^] consists to add in the configuration space two Dirac's 
functions to the solution. The coefficients of these 5 function are computed in 
order to satisfy the BC. The source S^ +1 ^ 2 {xk) is replaced by 

S{x k ) = S{x k y +l ' 2 + hS(x k , -1) + b 2 S(x k , 1) . (34) 

These three methods can be used either with explicit or implicit temporal 
scheme. From the numerical point of view, the implicit scheme must be treated 
carefully, especially for second order equations when the time step is large. 

4-5 Details about the boundary condition implementation 

Let us consider the heat equation (26) with v(x) = and fi(x) = 1 and let 
us suppose that the time discretization is performed according to a first order 
fully implicit time scheme. For At = 0.1, the system to be solved is then (cf. 
Eq. (16)): 
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where the R.H.S. of the system, S n , is the Lanczos modified source defined by 
Eq. (33) or by Eq. (34) depending on the choice of the BC treatment method. 
The simplest way to find the value of the coefficients b\ and b 2 is to solve first 
the system having as source the vector S n only (b\ = b 2 = 0). In this way one 
particular solution G)P ar of the system is found. The second step consists in 
finding two homogeneous solutions O^ 1 and 0^ 2 by solving twice the system 
(35) with the source 5n-i^ and o~N-2,n (for the Lanczos method) and finally 

5 In the book by Gottlieb & Orszag [21], this method is called collocation approx- 
imation (page 14 of [21]). We prefer the term dual Lanczos method for obvious 
reasons. 
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by computing b\ and 62 in order to fulfil the required BC. The same procedure 
holds for the dual Lanczos method. We recall that the Chebyshev coefficients of 
the Dirac functions at the interval boundaries, namely 5(x + 1) and 5(x — 1), 
are respectively 5 n = (— l)™ -1 and 5 n = 1. Of course, as already said, the 
matrix of the above system can be reduced with simple linear combinations to 
a penta-diagonal matrix and consequently finding solutions of a banded system 
having 3 different R.H.S. is a very little time consuming. The homogeneous 
solution can be used to satisfy the dual BC by imposing, for example, that 
the last coefficient of the solution vanishes. By looking at the matrix of the 
system (35), it appears that the extra-diagonal coefficients for At = .1 are 
much larger than the coefficients of the main diagonal. It follows that each 
solution of the system depends strongly on the value of the highest order 
coefficients of the source (this is a consequence of the fact that the problem 
is not well posed without BC). Moreover, small values of the highest order 
coefficients lead to very large solutions. When the linear combination of the 3 
solutions required to fulfil the correct BC is performed compensations occur 
to give the correct solution. These compensations are source of important 
round-off errors. In practice, At oc 1/iV 3 is a necessary condition to have high 
precision results. Another procedure works in the opposite case (At high) and 
is easy to implement for the Lanczos method. It is the following one: look for 
a solution of the first equation having some condition in the coefficient space, 
for example having the two first coefficients vanishing. This can be obtained 
easily by replacing the system (35) by 
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(36) 



In this way we have found a particular solution having the two first coefficients 
vanishing. Find then two other homogeneous solutions having as source 5^ n 
and #2,n and make a linear combination in order to fulfil the BC. It is obvious 
that the solution obtained in this way is the same than that obtained by 
solving the original system (35): both solutions satisfy indeed the first N — 2 
equations and the BC; the unicity theorem garanties that the two solutions 
are identical. 

The matrix of the new system has the coefficients of main diagonal which 
increase when N increases and the solution is computed with high accuracy 
for At > 1/iV 4 . In fact, in the opposite case, the determinant of the new 
matrix vanishes for At = and the round-off errors become important for too 
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small values of At. Therefore the solution method depends on the value of TV 
and At. We recommend the Lanczos method for large values of At and the 
dual method in the opposite case. For first order equations, or for the advective 
diffusive equation with high Peclet numbers, the dual Lanczos method is more 
convenient. 



4-6 Artificial spectral diffusivity 



In classical hydrodynamics with finite difference methods, artificial viscosity 
or diffusivity is often used in order to avoid singularities in the numerical 
solutions and/or to stabilize the numerical scheme. 

Within Fourier expansion, a high order operator is added to the RHS of the 
PDE in order to obtain a smoother solution. In the Fourier space, this artificial 
diffusivity is of the kind oc k 2L [2] (L integer number) where k is the spatial fre- 
quency. For L — 1, this term represents the ordinary diffusivity. For L — > 00, 
its tends to the evanescent viscosity introduced by Maday et al. [36]. Within 
Chebyshev expansion, spectral viscosity can be implemented via a power of 



[37]. This is an ex- 
for Fourier SM. Its 



the degenerated operator V = vT~— x 2 d/alx \/l — x 2 d/dx 
tension to Chebyshev SM of the technique developed in [2 
matrix representation is a diagonal matrix: V nm = —n 2 5 nm . Analogous linear 
artificial viscosity methods were used in numerical magneto-hydrodynamics 
for handling sharp gradients induced by non-linearities [29,52,35]. Non-linear 
artificial viscosity (hyperviscosity) has been introduced in [48] and applied to 
MHD problems [47]. 



5 Examples of time evolution 



5.1 Advection equation 



In view to compare the different approximations we start with the advec- 
tion equation (26) with fi(x) = and v(x) = 1. We choose Q(x, 0) = 
and 0(— l,t) = sin(47rt) 2 . The analytical solution of this PDE is Q aunSu {x,t) = 
sin(47r(t— x)) 2 and will be used to compute the numerical error. Spectral meth- 
ods are not well suited to solve hyperbolic equations; this example is worth 
to show why. Figure 3 shows the solution computed by means of a Chebyshev 
expansion and an implicit Crank-Nicolson scheme, with the Lanczos approxi- 
mation (Sect. 4.4). The corresponding numerical error is shown in Fig. 4. The 
solution is C 1 . Therefore the errors decreases as 1/N 2 and Gibbs phenomenon 
is heavily present. The solution becomes C°° only when the front of the wave 
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N = 65 




-1 -0.5 0.5 1 



Fig. 3. Profiles at different time steps (long dashed, short dashed and solid line 
successively) of the solution of the advection equation (Eq. (26) with v(x) = 1 and 
H = 0). The initial conditions are @(x, 0) = and the BC is 6(— 1, t) = sin(47rt) 2 . A 
sinusoidal wave propagates from the left to the right of the interval. The collocation 
points are shown on the plots (N = 65, At = 5 x 10~ 3 ). 

N = 65 

0.005 | ■ , ■ , ■ , ■ 1 




Fig. 4. Profiles at different time steps of the error of the numerical solution shown 
in Fig. 3. A Gibbs phenomenon appears when the wave is not arrived at the end of 
interval yet (see text). 

has crossed the right boundary. In the present case, A^ = 65 and At = 5x 10~ 3 . 
When the solution is C°°, the error is due to the temporal scheme which is 
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Fig. 5. Profile at different time steps of the Heaviside wave propagating from the 
left to the right and computed by means of Lanczos approximation. The solid line 
is the numerical solution and the dashed line the exact one. 

a second order implicit one. Decreasing At by a factor 2 reduces the error 
by a factor 4. From Fig. 4, it appears that the error increases linearly as 
the wave propagates. This is due to the implicit scheme which introduces a 
dissipative term, which gives an error on the amplitude of the wave. An explicit 
leap-frog scheme would avoid this drawback. However, whereas it is possible 
to use a leapfrog temporal scheme for Fourier expansion, such a scheme is 
unconditionally unstable for Chebyshev expansion (see [21] p. 109). Note that 
with N — 21, the error due to the computation of the spatial derivative is still 
less than the one due to the temporal scheme when the solution becomes C°° . 

We conclude that implicit methods should be avoided for solving hyperbolic 
equations when the phase must be computed with a high accuracy. Finally, 
it is interesting to compare the above results with the one obtained with 
finite difference methods. Such a comparison can be found in [21] p. 135. 
Unfortunately, these authors compare the results obtained with SM versus 
that obtained by finite difference methods by means of the £ 2 error and, 
therefore, do not give the error on the phase and amplitude of the wave. 

In the second example, we compare the solutions obtained with the the Lanc- 
zos and dual Lanczos approximations and study the effect of the spectral 
viscosity. Let us consider again the advection equation with same initial con- 
dition Q(x, 0) = but with a sudden excitation at x — — 1. The solution is a 
Heaviside function traveling from the left to the right. 

Figure 5 shows the solution computed with the Lanczos approximation. Note 
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Fig. 6. Same as in Fig. 5 but with the dual Lanczos approximation. 

that, in spite of a strong Gibbs phenomenon, the propagation velocity is cor- 
rectly computed. Figure 6 shows the solution computed with the dual Lanczos 
approximation (Sect. 4.4). The Gibbs phenomenon is much less important. 
Figure 7 shows the effect of the spectral viscosity (Sect. 4.6) (we have used 
a viscosity operator equivalent to V 2 ). The apodisation effect of the spectral 
viscosity is well visible. Note that, in spite of the spectral viscosity, the velocity 
of the front wave is still accurately computed. 

5.2 Boundary layers and non local boundary conditions 

In this section we give some result obtained by solving the equation (26) in 
order to show how a boundary layer can be described by means of SM and to 
show how non local BC works. Fig. 8 shows the solution at different time steps 
of the equation (26) with v(x) = 1 and /i(x) = .01. The initial conditions are 
the same as in the previous example: Q(x, 0) = 0. We suddenly switch the left 
BC to the value 1. The other BC is 0(1, t) = 0. The Lanczos approximation 
is used. The different plots of Fig. 8 show the wave propagating from the 
left to the right. The solution is similar to the one showed for the advection 
equation but, because of the diffusivity, the wave front is less stiff than in 
the case of pure advection. When the wave reaches the boundary, the layer 
is formed. Its thickness is of the order of 2/N p = 10™ 2 . Figure 9 shows the 
comparison between the solution fulfilling the above BC (solid line) and that 
with a non-local BC (dashed line). The coordinate x was rescaled in order to 
show the structure of the boundary layer. The BC are 0(— l,t) = 1 and the 
second one was chosen in order to minimize the sum of the square of the last 
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Fig. 7. Same as in Fig. 5 but with dual Lanczos approximation and spectral viscosity. 
The two profiles, beside the analytical solution (dashed line), correspond to two 
different values of the spectral viscosity. Note that the discontinuity is smeared out. 
In spite of the strong spectral diffusivity, the velocity of propagation of the shock is 
still correctly computed. 

three coefficients of the solution. Note that the new solution differs only inside 
the boundary layer (0.97 < x < 1). The difference between the two solutions 
is ~ 1 x ICT 12 for x < 0.5. With a non local BC, the numerical stability is 
achieved for any value of \i > 0. 



5.3 Kompaneets equation and space compactification 



The last example is the resolution of the Kompaneets equation ([41]). This 
equation describes the evolution of the photon distribution function in a bath 
of plasma at thermal equilibrium at temperature T within the Fokker-Planck 
approximation (see e.g. [41]). The Kompaneets equation reads: 



dAf _ d 
dt dv 



v 2 ^ - 2(3vAT + v 2 M + A/" 2 
av 



(37) 



where M{y, t) is the photon density (function of the frequency v (0 < v < oo) 
and the time t) and f3 is the plasma temperature in adimensional units. Its 
asymptotical solution is the Bose distribution: 



v 2 



N{y, oo) = -—3 r — - , (38) 

v ; exp(z///3 + //)-!' y 1 
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Peclet number = 200 

N=65 




Fig. 8. Profiles at different times of the solution of Eq. (26). Note the formation 
of the boundary layer near x = +1 (N=65, v(x) = 1, fi(x) = .01, N p = 200, 
At = .015). 
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Fig. 9. Same as in Fig. 8 with different boundary conditions. The full line is the 
solution of with two local B.C. and the dashed one is the solution with a non-local 
B.C. The x axis has been scaled in order to show the structure of the boundary 
layer. 

where \x > is the chemical potential. It is obvious that the Kompaneets 
equation preserves the number of the photons. If the photon density becomes 
larger than a critical value, then the photons condense at the frequency v — 
(Bose condensation). 
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To solve the Kompaneets equation, it turns out to be convenient to compactify 
the frequency space. This compactification can be performed by means of a 
new variable x = (u—1)/ (x E [—1,1], u — (l+x)/(l—x)). Equation (37) 
then becomes 



+ 



v 



dM 4(3v 2 d 2 N 
~dt ~ (1 + z/) 4 dx 2 ' (1 + z/) 2 

+ 2{v-f3)N 



+ 2M- 



2u 2 (3' 
1 + u 



dx 



(39) 



If JV(-1,0) = A/"(1,0) = 0, then N{-l,t) = Af(l,t) = for all t > and 
the coefficients of the derivatives vanish at the boundaries x = — 1 and x — 1. 
Therefore, no BC are required. Moreover, the Courant conditions are fulfilled 
for At oc 1/N 2 and implicitation is not necessary. 



6 3-D decomposition in spherical coordinates 



Spectral methods are well suited to handle the pseudo-singularities tied to 
spherical-like coordinates (r,9,(f>) at r = and sin 9 = 0. The various scalar, 
vector and tensor fields involved in any physical problem cannot be any arbi- 
trary function of these coordinates. Any quantity has to satisfy the so-called 
regularity conditions on the axis sin9 = and at the origin r = 0. These regu- 
larity conditions can be taken into account if an appropriated basis of spectral 
expansion is chosen. Even if regularity conditions cannot be called boundary 
conditions, this method to handle pseudo-singularities can be compared to the 
Galerkin approximation. 

As an example, let us consider a scalar function f(t,r,8,<f)), with r e [0,1], 
6 G [0, 7r], and <p G [0,27r[ where (r,9,(f>) are the usual spherical coordinates. 
Assuming that / is a regular scalar function means, in the context of the 
spectral approximation, / 6 C u , that is : 

f = J2a ijk x i y j z k , (40) 

where 

x = r sin 9 cos , y = r sin 6* sin and ^ = rcos^. (41) 

are the usual Cartesian coordinates. Now, substituting (41) into (40) shows 
that /, considered as a function of (r, 9, 0), reads 

f(r, 9,<t>)= a 'jk m r m+2j+k sm m+2j 9 cos k 9 exp(im0) . (42) 

j,k,m 



28 



In other words, this means that / is (obviously) a periodic function in the 
^-direction: 

fir, 6,<f)) = Y^ am(r, 9) exp(im0) , (43) 

m 

where the coefficients a m (r,9) satisfy 

a m (r, 9) = sin m 9 ]T a, m (r)P, m (cos 9) , (44) 
i 

where P™ are the associated Legendre functions. The coefficients a; m (r) satisfy 

ai m (r) =r l Y,ajimr 2j ■ (45) 
j 

A way to handle the regularity conditions on the axis sin6> = is to expand 
the angular part of / in a series of spherical harmonics 

f(r,9,<f ) ) = J2ai m (r)Y l m (9, ( f ) ) . (46) 

l,m 

Such an expansion ensures that any pseudo-singularity involving sin^ = 
is automatically handled. Moreover, spherical harmonics are eigenvectors of 
the Laplacian operator A. This property allows easy inversion of this opera- 
tor which appears in a lot of physical equations (particularly D'Alembertian 
operator, dissipative terms of hydrodynamics equations, Poisson equation for 
the gravitational field, see Sect. 7). 

However, from a numerical point of view, this expansion is not efficient in the 
sense that the transformation in the ^-direction is an expansion in associated 
Legendre functions P/ m (cos6>) which requires oc iV 2 operations, N being the 
number of degrees of freedom, compared to iV log N operations in the case of 
a Fourier transformation (thanks to the FFT algorithm). 

Moreover, it is to be noticed that the regularity conditions described above, 
namely that the coefficients a m (r,9) of the expansion (43) have to vanish on 
the axis as sin" 1 9, are not necessary to handle the singularities present in the 
differential operators. More precisely, it is easy to show that any singularity 
are automatically handled if the coefficients a m (r,9) satisfy 

a m (r,9) = sm M ^ m >°\9)Y,ai m (r)cos(l9) , (47) 

i 

where o is the order of the differential operator in the ^-direction. The same 
considerations hold for the expansion in the r-direction. 
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Fig. 10. Expansion of a 3-D scalar function in spherical-like coordinates (r,9,(p). 

For these reasons, we expand / as a truncated Fourier series in the ^-direction 
which requires oc log operations. We perform an analytical continuation 
of the coefficients a m (r,9) on the domain 6 G [— n, +n] and expand them in 
an even (resp. odd) Fourier series in the ^-direction for even (resp. odd) m. 
This transformation requires Ng log Ng operations. The resulting coefficients 
ct>im(r) are analytically continued on the domain r e [—1, +1] and expanded in 
a series of even (respectively odd) Chebyshev polynomials depending on the 
parity of /. This later transformation requires N r \ogN r operations. The above 
procedure is detailed in Fig. 10. 



7 Scalar Poisson equation in spherical-like coordinates 



We will describe methods which are specifically designed to solve PDE which 
contain scalar and vectorial differential operators in a space domain diffeomor- 
phic to a sphere and in a compactified space. All the numerical algorithms, but 
one, are fast algorithms. The only exception is the Legendre transformation 
which requires oc N r log N r NgN^, log instead of oc N r log N r N e log NgN^ log 
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operations. 



7. 1 Scalar Poisson Equation in a sphere 



Consider the Poisson-like equation in spherical coordinates, 

<9 2 $ 2<9$ 1 /<9 2 $ cosfl<9$ 1 <9 2 $\ _ c ff $ ^ , 48 n 



Solving by iteration this equation by means of spectral methods in a spherical 
region is described in [10]. Suppose that the source is known at the j th iter- 
ation step. The solution of the Poisson equation at the (j + l) th step can be 
obtained by solving the following system of ordinary equations on the spherical 
harmonics coefficients of S, say Si m (r): 

( d 2 2d + . . a . . 

{^ + - Tr - ± ^ 1 )^(r) = S lm (r), (49) 



where $; m are the spherical harmonics coefficients of The involved Legendre 
expansion is computed from the Chebyshev expansion in 9 by means of a 
matrix multiplication, as detailed in Appendix A.l of ref. [16]. 

Galerkin approximation is used to handle the singularity of Eq.(49) at r = 0. 
For / = the solution is expanded on even Chebyshev polynomials T 2n (r), 
(0 < r < 1, < n < N). For / even, the solution is expanded on the set 
of functions T 2n (r) — T 2n+2 (r). These functions as their first derivative vanish 
at r = 0. An analogous treatment is performed for odd values of I [16]. The 
matrix representation of the above operator can be easily reduced to a penta- 
diagonal one. This matrix has a vanishing determinant due to the existence 
of homogeneous solution $f m (r) = r l . 

The BC are satisfied by adding to a particular solution a homogeneous solution 
which must be computed in the coefficient space (see Sect. 4.5). Note 
that when the source does not depend on <p (2-D problems) there exists a 
fast algorithm for solving the Poisson equation that does not need a Legendre 
expansion (see Appendix A. 2 of [16]). The technical procedure outlined is 
here is described in details in Sect. 4.5, since the problem is reduced to a 
one-dimensional one for each value of (l,m). 
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7.2 Scalar Poisson Equation in a shell 



The solution of the Poisson equation in a shell R in < r < R out is performed 
in a very similar way. The main differences are: there is no singularity in 
the coefficients of the radial part of Eq. (49) and there are two homogeneous 
solutions &frn( r ) = rl an d = l/ r?+1 - The last analytical homogeneous 
solution can be used to fulfil the two BC. However, it turns out that better 
results are obtained by using an homogeneous solution obtained by means of 
the Lanczos approximation. In practice an aspect ratio R out /Ri n larger than 
2 should be avoided. 



7. 3 Solution of the Poisson equation in a compactified space 

Some of the Einstein equations lead to Poisson equations for which the source 
S(r, 9, 0) fills all the space and vanishes at least as oc 1/r 4 when r — > 00. Noting 
that the solution can be expanded in power series of the variable u — 1/r, it 
is natural to solve the Poisson equation in two domains D x and D 2 defined 
respectively by < r < 1 and 1 > u > and to match the two solutions 
and their derivatives at the common boundary. The radial part of the Poisson 
operator (with respect to u) in the domain D 2 reads: 



The source must thus be divided by the factor u 4 . This division must be 
performed in the coefficient space in order to avoid the singularity at u — 
and to minimize the roundoff errors. The reader will find all details in [12]. 
Extension of the above technique for more complicate cases (few domains) are 
straightforward. Generalization to parabolic equations is also straightforward. 
On the contrary, the above space compactification cannot be used for the wave 
equation. 



7.4 Poisson equation in a non-spherical domain 

Rotating or/and binary stars are not spherical but their surface are diffeomor- 
phic to a sphere. In order to fulfil the BC and to reduce truncation errors, it 
is useful to introduce a system of coordinates adapted to the geometry of the 
star. It is easy to find such a coordinates system, at least for star-like domains 
(cf. [12]). Let r = F(9, 4>) be the equation of the surface of the star. Consider 




(50) 
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the new coordinate system (£, 9', <p') defined by [12] 

r = fl( + e 3 F(^)+^i), 0' = 9, 0' = 0, 



(51) 



where F and F are such that (i) F(9, <p) = a + -F(#, 0) + F(9, <p) and (ii) the 
Fourier expansion of F(9,<f>) (resp. F(9,<f>)) with respect to contains only 
even (resp. odd) harmonics. The property (i) ensures that the surface of the 
star is given by £ = 1, whereas property (ii) is necessary to have a regular 
coordinates transformation. Under the coordinates transformation (51), the 
Poisson equation can be written in the following way: 

A A ' $ = Blk ++cl i^ + ' $j ' (52) 

where A' is the ordinary Laplacian with respect to the coordinates (£, , <j>') 
(i.e. A' has the same form as the operator (48) when (r, 9, 0) are replaced 
by (£,0'4>')) and A, B and C are functions of the coordinates (£,,9'<f)'). Equa- 
tion (52) is solved by means of an iterative procedure. Since A is not a con- 
stant, a procedure similar to that described in Sect. 4.1 (Eq. (29)) is used. 
No convergence problems were found. Note that the above method can be 
inefficient if the source S does not depend on $ and on its derivatives because 
the iterative procedure would still be required, whereas the equation would 
be linear in $. However, in the cases of interest, the source S is a non linear 
function of $ and of its first derivatives. Therefore an iterative procedure is 
required anyway. The computing time for solving Eq. (52) is similar to the 
one for solving Eq. (48). 



8 Vectorial equations in spherical-like coordinates 



The divergence and curl operators acting on a 3-vector V read: 



8Vr 



dVe cos 9 



V -V = — + - \ 2V r + — - + 

dr r \ 89 sin 9 



V e + 



1 dVs 



sm9 dcj) 



(53) 



and 



VAV 

Vav 

VAV 



r = [de(V <p sm9)-d^V e }/(rsm9) 
= [d <p V r -d r (rsm9V <p )}/(rsm9) 

u 

= [dr(rVe)-d e V r } /r , 



(54) 
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— * 

where V r , Vg, denote the components of a V with respect to the orthonor- 
mal basis (d r , 1/r dg, l/(rsin#) d^) associated with spherical coordinates 

It is worth to note a typical behaviour of the spherical components of a vector. 
Consider a constant vector V whose Cartesian components are V x — 1, V y — 
0, V z — . The corresponding spherical components are V r = sin 9 cos <fi, Vg = 
cos 9 cos 0, = — sin <f>. From the expression of the divergence (53), it is easy 
to see that there are three singular terms at r = and 9 = 0, n while the sum 
of these terms is regular. This is due to the fact that the spherical components 
of a well behaved vector are not independent. 

A way to overcome this problem is to compute and to add the terms which 
generate the singularity before dividing by r and sin# and to perform the 
division on the sum. Performing these operations in the coefficient space al- 
lows to resolve the indeterminations at the singular points and to reduce the 
roundoff errors. Another method consists in subtracting to V the constant 
component whose divergence vanishes. This decomposition can be easily done 
in the coefficient space. 

8.1 Vector decomposition 

The Clebsch theorem states that a vector V can be uniquely decomposed 
into a sum of a divergence-free vector W and a curl-free vector V $. This 
decomposition results from the resolution of a Poisson equation: the divergence 

— * — * — * — * — * 

of V = W + V $ gives indeed A$ = V • V, which has to be solved to get $. 

— * — * — * —* 

W is then computed according to W = V — V$. The above primitives are 
computed according to the method exposed in Sect. 3.2. Once a particular 
solution W p is obtained, a more general solution satisfying the desired BC can 
be obtained by adding to W p an homogeneous solution Wh = V $ha were <3>ha 
is an harmonic function. 

8.2 Solution of the equation V A V = S 

The curl operator VA being degenerated, the vectorial equation V A V = S 
can be solved only if the divergence of S vanishes (integrability condition). 

— * 

Because of this degeneracy, we can seek for a particular solution V for which 
one component vanishes. 

Taking for example V r = 0, the second equation of the system (54) V A V = S 
becomes d r (r V^) = —r Sg from which we get V$ — — £ J r r'Sgdr' . In the same 
way, we get Vg — - J r r'S^dr' . Once this particular solution is obtained, we 
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can add to it the gradient of some scalar function in order to obtain a more 
general solution satisfying the desired BC. 



8.3 Vectorial Poisson Equation 



The Einstein equations lead to a vectorial Poisson equation for the shift vector 
(Eq. (69) below): 

AV + XV(V-V) = S , (55) 



where A is some constant. Noting that A = —V A V A +VV- , Eq. (55) can 
be written as 

-V A V AV + (X + 1)V(V -V) = S . (56) 



We shall show how to find the solution of (56) in the spherical domain (0 < 
r < 1). Let us first consider the degenerated case A = — 1. The integrability 
condition is V • S — 0. Let us introduce the vector P = V A V which obeys 
to V A P = S. Following the method described in Sect. 8.2, we can solve this 
last equation and get P = P p + V(\I/ + ^ha) where P p is a particular solution, 
\& an arbitrary function and an harmonic function which can be used to 

— * — * 

fulfil the BC. The integrability condition imposes that A\& = —V • P p . The 
final equation for V is V A V = P + V(^ + *ha)- A particular solution V p of 
this equation (for \I/ha — 0) can be obtained by means of the method described 
in Sect. 8.2. The vector V g = V p + VS, where S is some arbitrary function, is 
again a solution of (55) but is not the most general one. Actually, we can add 
to it any solution of V A \4a = V^ha- Such a solution can easily be obtained 
by means of a spherical harmonics expansion of ^ha- 

The case A ^ -1 can be reduced to the previous one. Let us introduce V = 
W + V$ and apply the vector decomposition as described in section 8.1 to 
the source. $ can be computed by means of the divergence of both sides of 
Eq. (55). Now, the divergence- free counterpart W can be computed as in the 
previous case (A = —1). 



35 



8.4 Telegraph vectorial equation 



The methods described above to solve Poisson equations cannot be generalized 
to the parabolic or hyperbolic case. Consider the equation (telegraph equation) 

a ^¥ + + A ^ ' V) + a9 = § ' (57) 



where a, f3 and A are arbitrary constants. After decomposition of the source 
into a curl- free and a divergence-free component, the problem can be split in 
two parts. This is done mainly in solving a scalar telegraph equation for the 
curl-free component and a vectorial (with A = 0) telegraph equation for the 
divergence-free part. 

Eq. (57) explicitely reads: 



d 2 V r AdV r 1 / „d*Vr cos 9 dV r 1 d 2 V r 



D t y r _|_ " • r _|_ Zz_l. _)_ _ 2V r + + 

r dr 2 r dr r 2 \ r 89 2 sin9 89 sin9 2 dcj) 2 J 

--V.V = S r (58) 
r 

nT/ d 2 V e 2dV e 1 (d 2 V e , cos0 c>V 6 





-2^1^-^] =S e (59) 




Q r 2 r dr r 2 \ d9 2 sin 9 89 sin 9 

1 I ndV* dV r 2co80dV e 
r 2 sm6 \ 69 d(p sm9 op sm.9; 

where we have introduced D t V d = ad 2 V/dt 2 + (3dV /dt. It easy to see that 
these equations are badly coupled. Moreover, there is a lot of singular terms 
(cf. Sect. 8). The first equation is coupled with the others by V • V. Conse- 
quently, seeking for a divergence-free solution to the equation for V r reduces 
to a telegraph equation and can be solved mutatis mutandis by means of the 
method described in [5]. Once V r is known, we can seek for a solution of the 
system of coupled equation governing Vq, Va whose sources are 



Se = Se-2/r 2 d e V r 

S <j) = S <t> -2/(rsin9) 2 d <t> V r 



(61) 
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Introducing the couples of potentials (A, T) and (U, W), defined by 



OA 



1 or 



sin6> 



6., !('» + |)4 ( (62, 
r Vsm 9 acp 09 J 



1 fdU 

r 



1 dW" 



09 sin 9 dcj) 



V e , 



1 3U dW\ . , 

+ — I = V, (63) 



r \sin# 90 



allows to decouple Eqs. (59) and (60). 

Assume that the functions A and T are known by means of the resolution of 
the system (62,63). After some algebra, a solution of the telegraph equation 
(57) can be obtained by solving the associated equations 



D t U + AU = A 

D t W + AW = T , (64) 



where 



A = — 1 ( ^ I C0S6 9 J——) (65) 
dr 2 r 2 \d9 2 sin 9 09 sin# 2 dcj) 2 J 



It remains now to find the sources A and T of the system (64). Taking the 
angular divergence of S$ and allows to transform the system (61) as 

dS cos 9 - 1 dS^ 
A e ^A = — + — — -S e + . J , 66 
a9 sin 6* smfc* 2 ocp 



where Aq^ is the angular part of the ordinary Laplacian. This last equation 
can be easily solved by means of a Fourier-Legendre transformation. Once A 
is known, T can be easily computed. 

It is to be noticed that the use of Cartesian components V x , V y , V z would 
greatly simplify the resolution of the telegraph equation. The drawback of this 
short cut is that only some simple classes of BC can be easily implemented 
(e.g. V r = V e = 1/0 = 0). For more general BC, spherical components of V 

— * 

and S must be used. Finally, let us mention that the above method can be 
used to solve the vectorial Poisson equation too, since it is a special case of 
Eq. (57), obtained by setting a, (3 and A to zero. This integration method was 
used to solve a vectorial Poisson equation in studying the bifurcation points 
of general relativistic rotating neutron stars [5] (cf. Sect. 9.9 below). 
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9 Astrophysical applications 



9.1 Einstein equations 



In this section we give for completeness a short outline of the Einstein equa- 
tions [53], [58], [59]. This section is not strictly necessary to understand the 
following ones, and the reader that does not have some background of differ- 
ential geometry can skip it. 

Einstein equations describe the evolution of the metric g a p of a 4-D Rieman- 
nian manifold (Greek indices a, (5 run from to 3, Latin indices i,j run from 
1 to 3, x° = ct). It is convenient to write the distance element ds 2 on the 
following form: 

ds 2 = -a 2 {dx ) 2 + ^{dx 1 - (3 l dt) (dx j - (3 j dt) . (67) 



The 10 Einstein equations (1) are not independent because of the Bianchi 
identities (2). Therefore there is the possibility to impose 4 gauge conditions. 
If the maximum slicing and minimal distortion gauge is chosen (cf. [53]), the 
Einstein equations (1) can be written in the following symbolic way (where 
we have introduce the Cartesian components of 7^) 

Aa = p + Q (68) 

AA + =Si + Qi (69) 

Alog 7 = P + Q (70) 
1 <9 2 7 ifc 



c 2 dt 2 



&lik = Pik + Qik , (71) 



where 7 is the determinant of the space metric 7^, p, Si, P and Pik are re- 
lated to the matter energy momentum tensor T a p, A is the ordinary flat-space 
Laplacian and Q , Qi, Q and Q ik represent a formal source containing all 
the non-linear terms of the Einstein equations. It is easy to see that the first 
five equations are elliptical, the other ones are hyperbolic. The above equa- 
tions can simplify enormously if some symmetries are present. For example, in 
spherical symmetry, a non-static solution of the system (68)- (71) can be found 
by solving a system of two PDEs for the metric and two PDEs for the matter. 
An axial symmetric circular steady state solution can be found by solving a 
systems of 4 elliptical equations [14]. The problems becomes more complicated 
in the noncircular case. In this case, in spite of the axial symmetry, all the 
10 equations of the system must be solved [24]. In the next sections we shall 
describe some astrophysical applications obtained by our group by solving the 
Einstein equations (68)-(71). 
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Fig. 11. Time evolution of the stellar radius for four initial configurations close to 
the neutron star maximum mass [23]. The relative differences with the maximum 
mass are respectively -3.3 x 1(T 3 , -5.0 x 10" 5 , -3.6 x 10~ 6 and -3.8 x 1(T 6 
for configurations labelled A, B, C and D. The models A, B and C are on the 
stable branch of the mass-central density curve, whereas model D is on the unstable 
branch. Note that the high frequency small oscillations are real and corresponds to 
stable acoustic modes. They keep the same amplitude and same frequency for the 
four configurations. On the contrary, the real part of the frequency of the unstable 
Jeans mode tends to zero along the sequence A, B, C and its amplitude increases. 
For configuration D the frequency becomes imaginary, leading to a collapse (not 
shown on the figure for the benefit of presentation). 

9.2 1-D gravitational collapse of a neutron star towards a black hole 

Accreting neutron stars (e.g. in X-ray binary systems) can reach the maximum 
mass allowed by general relativity [46] , which marks the stability limit against 
radial perturbations [30] . Any subsequent accretion will trigger the collapse to- 
wards a black hole. This phenomenon may be an important source of neutrinos 
[25] . We have studied the collapse by means of a general relativistic and spheri- 
cally symmetric hydrodynamical code [22] , [23] . The initial configurations were 
solutions of the hydrostatic equations (Tolmann-Oppenheimer-Volkoff equa- 
tions). The high accuracy of spectral methods allows to take initial stable 
configurations, the mass of which differs from the critical mass by a few 10~ 6 
solar mass. The numerical round-off error cause the solution not to remain ex- 
actly static. It performs instead gentle oscillations of relative amplitude 1CT 10 
except for the fundamental mode, which grows exponentially for initial un- 
stable configurations (cf. Fig 11). The collapse is then computed up to the 
formation of an apparent horizon which leads to a coordinate singularity [22]. 
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9.3 1-D gravitational collapse within tensor-scalar gravitational theory 

We have recently studied the gravitational wave emission from 1-D gravita- 
tional collapse by solving the complete tensor-scalar and hydrodynamic equa- 
tions for a self-gravitating perfect fluid [43,44]. The initial conditions describe 
unstable-equilibrium neutron star configuration. In the case of shock forma- 
tion, we have coupled our Poisson code to a Godunov-type solver for the hydro 
given by Valencia's group [40], [31], [43]. We found that these kind of sources 
are not likely to be observed by future laser interferometric gravitational wave 
detectors (such as VIRGO or LIGO) if they are located at more than a few 
100 kpc. However, spontaneous scalarization could be constrained if such a 
gravitational collapse is detected by its quadrupolar gravitational signal since 
this latter is quite lower than the monopolar one. 

9.4 2-D Kerr black hole 

We have tested our method for solving the axisymmetric stationary Einstein 
equations against a non trivial 2-D analytical solution, namely the Kerr so- 
lution (see e.g. [56]), which describes rotating black holes. For treating this 
problem, four equations must be solved: Eq. (68) for a, only one of the three 
Eqs. (69) (that for (3^) and two equations equivalent to the system (70)-(71). 
These equations are solved outside the event horizon (given by the (isotropic) 
radial coordinate r = 1) up to infinity, by means of a compactification as ex- 
plained in Sect. 7.3 [6,5]. The horizon is a singular point for the coefficient of 
the operator in Eq. (69). This problem is a typical Von Neumann-Dirichlet 
problem. Figure 12 shows the error on the shift (3^ for a Kerr black hole with 
an angular momentum parameter a/M = 0.99. Note that for a/M = 1 the 
horizon shrinks to a point, so that the singularity becomes essential. 

9.5 2-D magnetized rotating neutron stars 

We have developed a numerical code for solving the coupled Einstein-Maxwell 
equations describing steady states of magnetized rotating neutron stars. From 
the mathematical point of the view, this problem amounts to solving four 
non-linear elliptical gravitational field equations, as in Sect. 9.4, one elliptic 
equation for the magnetic vector potential (the Grad-Shafranov equation) and 
one elliptic equation for the electric scalar potential. These latter two equa- 
tions are coupled to the gravitational field ones via the metric tensor. The 
steady-state configuration of the fluid is obtained by using a first integral of 
motion [3] . This code has been applied to compute the magnetic field induced 
distortion of rotating neutron stars (Figs. 13 and 14) in order to estimate the 
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Fig. 12. Shift vector component fift (top) and error with respect to the analyt- 
ical solution (bottom) for a Kerr black hole with angular momentum parameter 
a/M = 0.99 [5]. The different curves corresponds to different values of the 6 angle. 
The asterisks indicate sub-domain boundaries. 

resulting gravitational radiation [7]. Note that a direct comparison study has 
been conducted between the spectral code and two codes based on different 
methods (finite differences) in the case of non-magnetized rapidly rotating 
neutron stars [45]. 

9.6 2-D hot new born neutron star 

We have computed models of differentially rotating proto-neutron stars using 
realistic equations of state of dense hot matter in the framework of general 
relativity [27] [28]. We found a minimum period of uniform rotation of cooled 
neutron stars which formed directly (i.e. without a subsequent significant ac- 
cretion of mass) from proto-neutron stars with shocked envelope of about 

1.7 ms. This strengthens the hypothesis that millisecond pulsars are accretion 
accelerated neutron stars. 
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Fig. 13. Magnetic field lines and surface of the star (thick line) for a neutron star 
deformed by a huge magnetic pressure [3]. 



Magnetic field lines 
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Fig. 14. Magnetic field lines and surface of the star (thick line) for a realistic 
magnetized neutron star, with superconducting interior [7] . 

9. 7 3-D stellar core collapse 

A fully 3-D hydrodynamical Newtonian code for self-gravitating bodies has 
been developed within spectral methods and employed to compute the gravi- 
tational wave emission from the infall phase of a type II supernova [17]. The 
initial conditions are self-consistent configurations of rotating stellar cores em- 
bedded in the tidal field of a companion, hence are fully triaxial. They are 
computed by means of an iterative procedure. 
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Fig. 15. Shape of the grid (adapted to the star) at different times during a close 
encounter with a massive black hole (the grid is initially the fourth of a sphere) [39] . 
This computation has been performed with only 17 degrees of freedom in r, 9 in 9 
and 8 in cf>. 

9. 8 3-D tidal disruption of a star by a massive black hole 

We have studied the close encounter of a star and massive black hole, possibly 
at the centre of a galaxy [39]. Tidal forces compress dramatically the star, as 
shown in Figs. 15 and 16. This phenomenon has been proposed as a possible 
mechanism for gamma-ray burst generation [19]. The 3-D Euler equations 
are solved in the external field of the black hole. As already said in Sect. 1, 
this application demonstrates the capability of SM in computational resource 
saving. 



9.9 Spontaneous symmetry breaking of a rapidly rotating neutron star 

We have computed the location of the bifurcation point between (axisymme- 
tric) MacLaurin-like and (triaxial) Jacobi-like configurations along sequences 
of rapidly rotating stars in the framework of general relativity [4,5]. To tackle 
this problem, we solved Eq. (68), the three Eqs. (69) (by means of the tech- 
nique described in Sect. 8.4), and two equations like (69) and (71). The grav- 
itational radiation is neglected. The matter distribution is obtained thanks 
to a first integral of motion and depends on the equation of state of nuclear 
matter. Such mechanism can be a powerful source of continuous gravitational 
waves from accreting neutron stars in binary systems [8,9]. 
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Fig. 16. Time evolution of the central density of star during its trip around a 
massive black hole [39]. 




Fig. 17. Rapidly rotating neutron star, which has become triaxial due to the viscos- 
ity driven instability [4,5]. The ellipticity in the equatorial plane is about 0.1. Notice 
the cusp at the stellar equator in direction of the semi-major axis (x-axis), where 
the rotation is Keplerian, and which is absent along the semi-minor axis (y— axis). 

9.10 Close binary systems of compact objects 



We have developed a numerical code for computing quasi-equilibrium configu- 
rations of relativistic binary systems [12]. This code is based on a multi-domain 
spectral method, i.e. the whole space is divided into various domains and a 
spectral method is used in each domain, as described in Sect. 7. The external 
domain extends to infinity thanks to some compactification. The boundaries 
of each domain are chosen in order to coincide with physical limit surfaces, 
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Fig. 18. Logarithm of the relative global error of the numerical solution with respect 
to the number of degrees of freedom in 9 for a Roche ellipsoid for an equal mass 
binary system and Q 2 /(7rGp) = 0.1147 (the numbers of degrees of freedom in the 
other directions are N r = 2Ng — 1 and N v = Nq — 1) [12]. Also shown is the error 
in the verification of the virial theorem. 

such as the surface of a star (cf. Sect. 7.4). In this manner, the physical discon- 
tinuities in some fields or their derivatives are located at the boundary of the 
domains and the applied spectral methods are free of any Gibbs-phenomenon, 
resulting in a very high precision. We use a spherical-type coordinate system 
(r', 6', <$') which maps the physical domain to the unit sphere. This technique 
is able to treat any physical boundary which is starlike; this includes any re- 
alistic shape taken by the surface of a rotating neutron star in a tidal field. 
Among the tests passed by the code, the comparison with analytic solutions 
(MacLaurin and Roche ellipsoids) shows that the achieved precision is of the 
order 10~ 10 (cf. Figure 18). The final stages of NS binaries are expected to be 
irrotational with respect to an inertial frame (see [11] and references therein). 
This means that when seen in the co-orbiting frame, each star is counter- 
rotating with respect to the orbital motion. The velocity field with respect 
to the orbiting frame must then to be computed. A typical numerical solu- 
tion is shown in Fig. 19. The first numerical results about irrotational binary 
configurations in general relativity have been presented in Ref. [13]. 



10 Conclusion 

We hope to have convinced the reader that spectral methods are a highly 
valuable tool in the field of relativistic astrophysics. They can lead to results 
typically several orders of magnitude more accurate than those provided by 
finite difference methods. This high accuracy is very useful for studying sta- 
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Fig. 19. Velocity field with respect to the orbiting frame, inside a star which is 
member of a close binary system. The plane of the figure is the orbital plane. The 
companion is located at the right. 

bility problems (Sect. 9.2 and 9.9). Of course, spectral methods have their 
drawbacks. Let us summarize the advantages and drawbacks from our expe- 
rience: 

• Advantages: 

■ very high accuracy for treating smooth fields and their derivatives (evanes- 
cent error); 

• economy of number of degrees of freedom (grid points); 

• rigorous treatment of boundary conditions; 

• multi-domain (multi-grid) technique easily implemented without any loss 
of accuracy; 

• rigorous treatment of regularity conditions associated with singular coor- 
dinates such as spherical coordinates; 

• efficient algorithms for von Neumann-Dirichlet problems; 

• modularity; 

• well developed mathematical theory; 

• lack of robustness (this is an advantage because no solution can be found 
for ill-posed problems); 

• Drawbacks: 

• spurious oscillations when treating discontinuous fields (Gibbs phenomenon); 

• severe Courant conditions and therefore need of implicitation (Sect. 4.1); 
Note however that in the 2-D Fourier case, spectral hyperviscosity has 
proved to be efficient in shock handling [48,47]. 

Clearly there is no ideal numerical method for treating all the problems and 
a combination of different methods for solving a given problem may turn out 
to be fruitful, as preliminary results presented in Sect. 9.3 indicate. It would 
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also be desirable to compare results obtained by spectral methods with that 
obtained by other methods, as was done for the problem of rotating neutron 
stars (cf. Sect. 9.5) and extensively for hydrodynamical shocks [57]. 

In our opinion, one of the main advantage of spectral methods is that they 
allow a rigorous treatment of the regularity conditions associated with the 
singularities of spherical coordinates. This is important because spherical-like 
coordinates are obviously much more adapted than Cartesian coordinates for 
describing objects like stars or black holes. In particular, the number of points 
required to treat a star with Cartesian coordinates with a reasonable accuracy 
is considerably higher than that required by spherical coordinates. Moreover, 
the surface of the star (or the black hole horizon), where boundary conditions 
may have to be set, is described more simply with spherical coordinates, espe- 
cially when surface-fitted spherical coordinates are introduced as in Sect. 7.4. 
Also the boundary conditions for Poisson-type equations are better expressed 
in spherical coordinates than in Cartesian ones. 

We thank Prof. Jose-Maria Ibahez and an anonymous referee for their carefull 
reading of the manuscript. 
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